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Abstract 

The properties of the solvent affect the behavior of the solution. We propose a model that accounts 
for the contribution of the solvent free energy to the free energy of globular proteins in solution. For 
the case of an attractive square well potential, we obtain an exact mapping of the phase diagram of 
this model without solvent to the model that includes the solute-solvent contribution. In particular 
we find for appropriate choices of parameters upper critical points, lower critical points and even 
closed loops with both upper and lower critical points, similar to one found before In the general 
case of systems whose interactions are not attractive square wells, this mapping procedure can be 
a first approximation to understand the phase diagram in the presence of solvent. We also present 
simulation results for both the square well model and a modified Lennard- Jones model. 

1 Introduction 

In recent years there has been an enormous increase in the number of proteins that can be isolated, due 
to the rapid advances in biotechnology. However, the determination of the function of these proteins has 
been slowed by the difficulty of determining their crystal structure by standard X-ray crystallography. A 
major problem is that it is difficult to grow good quality protein crystals. Experiments have clearly shown 
that this crystallization depends sensitively on the physical factors of the initial solution of proteins. An 
important observation was made by George and Wilson [5], who showed that x-ray quality globular 
protein crystals only result when the second virial coefficient, B2, of the osmotic pressure of the protein 
in solution lies within a narrow range. This corresponds to a rather narrow temperature window. For 
large positive B2, crystallization does not occur on observable time scales, whereas for large negative B2, 
amorphous precipitation occurs. Rosenbaum, Zamora and Zukoski then showed |3] that crystallization of 
globular proteins could be explained as arising from attractive interactions whose range is small compared 
with the molecule's diameter (corresponding to the narrow window of -E^)- In this case the gas- fluid 
coexistence curve is in a metastable region below the liquidus-solidus coexistence lines, terminating in a 
metastable critical point. 

When a system undergoes a phase transition the change in the Gibbs free energy AG consists of 
two terms - an enthalpy change AH and and an entropy change AS*, with AG = AH — TAS. The 
change in the Gibbs free energy must be negative in order for the transition to occur. The enthalpy 
change is negative because in the separated state the more dense phase has a larger number of contacts 
and therefore its contact energy is lower. However, the dense state has a lower entropy and its entropy 
change is also negative. These two terms compete and at low temperatures the free energy change due 
to particles going from the dilute phase to the dense phase is negative and phase separation occurs. 
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However, above some temperature the entropy loss exceeds the enthalpy loss, so that the change in the 
Gibbs free energy is positive and phase separation does not occur. 

The solvent can change this picture dramatically 0] . Consider the free energy change of the particle 
going from the dilute phase to the dense phase. 

AG = AH so i ute — TAS so i u te + AH so i vent — TAS so i ven t (1) 

Here AH ao i ute is the enthalpy change of the particle between the dense and dilute phases, AS so i u t e is the 
entropy of the solute particle change, AH so i vent is the enthalpy change of the solvent and AS solvent is 
the entropy change of the solvent. The first two terms on the right side of (JIJ behave exactly as described 
above. The last two terms, however can be either negative or positive. 

As an example of this, consider water. Water has a large number of strong intermolecular hydrogen 
bonds in the bulk state. When one adds a hydrophobic particle to water, some of these bonds break. 
This leads to an increase in the entropy of water and to a decrease in its enthalpy. However, the 
further rearrangement of the water molecules around the solute particle leads to even stronger bonds, 
thereby decreasing solvent entropy and decreasing the solvent enthalpy. When particles aggregate the 
water molecules return to the original bulk state so the entropy and the enthalpy of the water increase 
to their original values |S1 H3 IEj< This is an example how the AS so ivent and the AH so i ven t terms in 

can be positive. If the sign of the total enthalpy and the total entropy change is positive, then 
decreasing temperature can lead to the opposite behavior, when the system tends to be separated above 
some temperature and uniform below. This leads to a lower critical point. This phenomenon has been 
experimentally observed in different systems 1101 1111 E] ■ The most interesting example of a lower 
critical point in context of our work is the phase diagram of sickle hemoglobin, which is thought to have 
an upside down fluid fluid coexistence curve |T3], as inferred from the observation of an upside down 
spinodal curve. 

Below the lower critical point the water molecules build strong hydrogen bonds around solute particles 
preventing them from aggregating. As one increases the temperature the free energy decrease of the 
aggregating particles surpasses the free energy gain of the water and the system phase separates. If one 
increases the temperature further, the total entropy loss due to aggregation becomes larger than the 
total enthalpy loss, so the phase separation becomes less favorable again and above some upper critical 
temperature the system becomes uniform. Thus for the hydrophobic particles one can have both an 
upper and a lower critical point; in other words, the fluid-fluid phase diagram has the form of a closed 
loop. This can be the case for different protein-water solutions. 

In this paper we present a simple model for the role of the solvent. We model the multicomponcnt 
protein-solvent system as a binary system in which the solute molecule is much bigger than the solvent. 
We assume that the role of other components in the system are subsumed in the effective solute-solute 
interaction. The role of a salt, for example, is assumed simply to screen the charge on the proteins, 
inducing the effective inter particle attraction, whose strength is controlled by the salt concentration. 
This assumption permits us to predict three qualitatively different effects of the role of the solvent for 
globular proteins in solution, all of which are consistent with experimental observation 0J^]. The first is 
a negative enthalpy and entropy of crystallization, which gives a normal liquid-solid line. Lysozyme is an 
example of this, as it has a negative enthalpy of crystallization, of the order of minus 75 k J . The second 
is a positive enthalpy of crystallization, which gives the upside down (retrograde) behavior of coexistence 
curves. This has been seen in hemoglobin C (CO-HbC), which has a high positive enthalpy of 155 fcj -1 
of crystallization. This system exhibits a strong retrograde solubility dependence on temperature (figure 
^) . Another example of retrograde behavior is chymotrypsinogen A JS] . The third is a zero enthalpy of 
crystallization, which gives vertical coexistence curves. Apoferritin has an enthalpy of crystallization close 
to zero and has a solubility which is independent of temperature (vertical behavior in the temperature- 
solubility diagram). 

We consider a situation in which the solute particle-particle interactions are described by a repulsive 
hard core together with an attractive interaction. For example, the interparticle potential might be 
described by the DLVO potential [TJ], which consists of a Debye-H-itckel interaction and a van der Waals 
interaction. Other model potentials are the square well and modified Lennard- Jones models that have 
been used successfully in recent years to describe globular proteins in solution. If the attractive part is 
short-ranged then these effective potentials can qualitatively describe the canonical phase diagram of the 
globular protein solution. We can therefore consider that solvent contribution was accounted in these 
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Figure 1: Solubility curve of CO HbC at conditions indicated at the plot. Points are the experimental 
results and the curve is fit using AH — 155k J /mol Details of the fitting procedure can be found in 0]. 
The figure was nicely granted by P. G. Vekilov. 

effective potentials and there is no need to add the solvent-solute interactions into potential. Examples 
of proteins that have canonical phase diagram are Lysozyme |17l I15| and 7D-crystallin |18j . However 
for proteins that have completely different phase diagram (like apoferritin or HbS, HbC) these effective 
potentials fail to describe the phase diagram even qualitatively. In this case the solvent contribution 
wasn't successfully included into effective potential. The model presented in this paper is an attempt to 
account for the solvent effect to describe the phase diagrams of the proteins like apoferritin, HbS, HbC 
at least qualitatively. 

We then include a simplified solute-solvent interaction to the solute-solute interaction, which is similar 
to one used recently to describe hydrophobic interactions . Our particular interest is globular proteins 
in solution, but our treatment in principle includes other systems. 

The outline of the paper is as follows. Section 2 contains a description of our model, while section 
3 contains a discussion of the particular case of a square well system with solvent. For this particular 
choice of solute-solute interaction, we are able to obtain the phase diagram of the square-well potential 
model with solvent contribution from the phase diagram of the square-well potential without solvent 
contribution. Section 4 discusses three cases of this mapping, for different choices of parameter values, 
that yield upper critical points, lower critical points, and closed loops, respectively. Section 5 contains 
a derivation of a modified Clausius-Clapeyron equation in the presence of the solvent, that takes into 
account the temperature dependence of the effective interaction potential. Section 6 contains the results 
of a simulation of a modified Lennard- Jones model plus solvent. The modified Lcnnard- Jones model 
was originally introduced to describe the phase diagram and nucleation rates for globular proteins in 
solution PHI ■ We show that for a particular choice of parameter values the solvent produces a fluid- fluid 
coexistence curve with a lower critical point. The liquid-solid coexistence lines for this model are also 
shown. Finally, section 7 contains a brief conclusion. 

2 Simple model for solvent-solute interaction 

In this section we will derive a phenomenological model that describes the influence of the solvent on 
a system of interacting protein particles. Consider a system of protein particles interacting via a short 
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ranged pairwise potential U(fij). The total energy of such a system is 

Uo(r If ) = ±Y, U °(\ T i- r i\) (2) 

where the subscript denotes the potential energy in the absence of the solvent. The corresponding 
Helmholtz free energy F (N, V, T) of the system is 

F (N,V,T) = -k B T\n (jWJfi J dr N exp(-/3[/ (r Ar )) (3) 

where N is the number of particles in a volume V, A is the thermal de Broglie wavelength, r the 
coordinates of the particles and k B is a Boltzman constant. To take into account the effect of the solvent 
on this system, we note that a particular solvent molecule can be either in the bulk or in a " shell" region 
around a protein particle (see Fig. |2J). We define the energy and entropy (in units of k B ) differences 
between these two states to be e w and As Wl respectively. In general these quantities depend on the 
temperature 11010]; i- e - e w — £w(T), As w — As w (T). Thus the free energy difference between a solvent 
molecule in the bulk and in the shell is Af = e w — k B TAs w . The total energy of the protein-water 
system therefore is 

TV 

U = U a + J2(^- k B TAs w )n$ (r N ) (4) 

where riw is the number of water molecules around the i th particle. We now assume that this number is 
proportional to the available area on the protein molecule surface and approximate this area as propor- 
tional to the maximum possible number of contacts of the protein molecule with other protein molecules. 
That is, the available energy is proportional to n c minus the actual number of contacts rip\ 

n$(r N )cxn c -nW(r N ) (5) 

We absorb the proportionality coefficient into e w and As w . Next, define two particles to be in contact 
when the distance between their centers of mass is less than some value X c a and greater than a, where 
a is the hard core diameter. Define the function 7(|i"j — r 7 - 1 ) to be equal to one if the distance between 
particles i and j is less than X c a, but greater than a and is zero otherwise. If the distance between the 
particles is less than their hard core diameter, the total energy of the system is infinite, which is the value 
of Uo(\ri — Tj\ < a), so one can define 7(|r.; — r,-| < a) = at that distance. The number of contacts for 
the i th particle is then given by 

n« = $>(l r *- r iD (6) 
Combining equations (@J , JSJ) and © , we find that the total energy of the system is 

= \ E ^(ki - r,-|) + (e w - k B TAs w ) I Nn c - ^J 7 ( r, r, ) | (7) 

We can combine the first and the last terms on the right hand side and define a modified intermolecular 
interaction parameter 

tf(|ri-r 3 -|) = l7o(|ri-r i |)-2(e to -A!BTAs to )7(|r i -r i |) (8) 

In summary, our model assumes that a water molecule can be either in the bulk or close to the protein 
particle (shell) ; the number of water molecules near the particle is proportional to the available area on 
the surface of the particle; the available area on the protein surface is proportional to the maximum 
possible number of protein-protein contacts minus the actual number of contacts. 
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3 Square well system with solvent 



3.1 Implementation of the general theory 

In order to simplify this expression for the energy, we consider a system with a repulsive hard core and 
an attractive square-well (SW) potential, whose range of attraction is A c . We can write Uq in the form 
of Uq = — 1/2 J^ili e o n p\ where e is the well depth. The total energy of such a system is 



1 N 

U{r N ) = --^(eo + 2e w - 2k B TAs w )n^ + (e w - k B TAs w )Nn c (9) 

i=l 

The first term on the right hand side is just the energy of the square-well potential with the well depth 
e = e + 2e w — 2k B T As w . This simplification allows us to map the phase diagram of this square- well 
potential model with solvent contribution onto the phase diagram of the square- well potential without 
solvent contribution. Denote the free energy of the square- well system with well depth e as Fo(N, V, T; e). 
The second term on the right hand side of I© does not depend on the position of the particles. Therefore 
the free energy of the SW system with solvent 

F(N, V, T) = F (N, V, T; e) + (e w - k B TAs w )Nn c (10) 

The free energy density of the system is 

f(p, T ) = foiP, T; e) + pn c (e w - k B TAs w ) (11) 

where /o is the free energy density of the square- well system (without solvent). Denote po and ujq as the 
chemical potential and grand canonical free energy density respectively. In this case pa = dfo/dp and 
= /o ~ PoP- To obtain the phase diagram for the square-well system one needs to solve the phase 
coexistence conditions po(pi,T) = po(p2,T) and loq{pi,T) — luq(p2 1 T). Suppose that we have obtained 
this phase diagram in some way (numerically, for example). The square- well coexistence curve has the 
following functional form 

- Tip) (12) 

e 

To obtain the phase diagram of a square well system with the solvent we should write the chemical 
potential and grand canonical free energy in terms of po and uiq. 

p{p, T) = p (p, T; e) + n c (e w - k B TAs w ) 

v(p,T) = u (p,T;e) (13) 

As we can see the phase diagram for the protein solvent system maps to the square well phase diagram 
with the well depth eo + — 2kTAs w : 

kT< r(p) (14) 



(^o 2c w 2kT coex As u 

Assuming that e w and As w do not depend on temperature, the relationship between the phase diagrams 
for these two models is given by 

k B T coex = rrT T (p) (15 

(1 + 2As w T(p)) 

In summary, we are able to obtain the phase diagram for an interacting protein system that includes 
the effect of solvent in terms of the interacting protein system without solvent, assuming that the protein- 
water attractive interactions are given by a square well potential; that the range of the protein-water 
interaction is equal to the range of the protein-protein interaction A = A c and that the parameters e w 
and As,,, are constants. 
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4 Upper and lower critical points and closed loops 



4.1 Upper and lower critical points 

In this section we consider examples of the phase diagram mapping for the square well model with solvent. 
Let r(p) in 112f) be a fluid-fluid coexistence curve. The critical point can be determined by the following 
condition 

RlcT 

=0 (16) 
dr{p)/dp 



$Pcoex 

If we define the left hand side of {H} as G(kT) then 

OkT coex 1 3 Pcoex — 



dG{kT)/dkT 



This means that the critical density of the modified system is the same as the critical density of the 
original square well system. The second derivative shows whether there is an upper or lower critical 
point. Taking into account that dr/dp = at the critical point, we obtain 

d 2 kT = d 2 r/dp 2 

dp 2 dG/dkT crit 1 > 

The original square well model has an upper critical point, so d 2 r(p cr ) / 'dp 2 < 0. Therefore the sign of 
the denominator determines the type of the critical point. If dG{kT cr ) / dkT > then there is an upper 
critical point. 

In our simplified model with constant e w and As w we can calculate dG/dkT explicitly dG/dkT = 
(eo + 2e w ) / (e 2 ). This gives the first condition for upper or lower critical points. The temperature in (|15|) 
should be positive. This gives the second condition required for having upper or lower critical points. 

e w > — eo/2, As w > eo/2kT® r - upper critical point (18) 
e w < —cq/2, As w < eo/2fcT c ° r - lower critical point (19) 

where T c ° r is the critical temperature of the original system without solvent. 
4.2 Closed loop phase diagram 

One of the possible fluid-fluid coexistence curves is a closed loop pp with both upper and lower critical 
points. In order to obtain this behavior from our model we must have temperature dependent e w and 
As w . We use the four level Muller, Lee, and Graziano (MLG) model of water [2111221 H] to obtain e w (kT) 
and As w (kT). In this model water molecules can be either in disordered (broken hydrogen bonds) or 
ordered states. When the protein particles are added, each state of a water molecule splits into shell and 
bulk states. Thus we have a four level model (Fig. [2J. The values of the energy levels and degeneracies 
used in [T] are following : E os = -2, E ob = -1, E db = 1, E ds = 1.8, q os = 1, q ob = 10, q d b = 40 and 
q ds = 49. Moelbert and De Los Rios considered a lattice model using the four level approximation 
and found by Monte Carlo simulation a fluid-fluid coexistence line in the form of a closed loop 

Using the above values of energies and degeneracies we can calculate the energy levels of the shell 
and bulk states: 

E -4- Ea P -P AE ° 

E < - J + e -p*E. ( 20 ) 

E * 1 + e -PAE b ( 21 ) 

In the same way we determine the values of the entropies for the bulk and shell states: 

-Lft-M fctS*^ ) (22, 
k~S b = ln( 1 + e _ 0AEb ) (23) 
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Figure 2: Energy levels for the MLG model for water. The configurations of the different states are 
shown schematically. The figure was nicely granted by S. Moelbert [25] . 



where AE S = Ed s — E os and AE b = — E ob . The parameters of our model e w and As w are the 
differences of the energies and entropies in the bulk and shell states respectively, 

e w (kT) = E s -E b (24) 
As w {kT) = S s - S b (25) 

Substituting (|24|l and l|25|l into l|14|l we can numerically map the square well phase diagram onto the 
MLJ-model phase diagram. Figure |21 shows the dependence of the left hand side of l|14|) versus the 
temperature of the system (mapping curve) in the presence of solvent for three different values of the 
protein-protein interaction strength. The horizontal line in fig. [3] represents the critical temperature 
of the original square- well system without solvent. The intersection of the mapping curve with the 
horizontal line is the critical temperature of the system in the presence of solvent. For each curve in 
fig. |3 there are two intersections with the horizontal line. Therefore the mapping procedure yields two 
critical points - one lower and one upper. The mapping of the entire original coexistence curve gives a 
closed loop type phase diagram for the MLG-type system. 

Using l|2(J|) - (|23|) we can obtain the low and high temperature values of e w and As w . This permits 
us to determine the upper and lower critical points. In the low temperature limit 

e w {0) =E os -E ob (26) 
As w (0) = ln(q os /q ob ) 



and in the high temperature limit: 



/ s E os + E ds E ob + E db 
<M°°) = 2 2 ^ ' 

As w (oo) = ln( ) 

qob + Qdb 

One necessary (but not sufficient) condition to have both upper and lower critical points is l|26|) to satisfy 
conditions (|I9|) and (|27|l to satisfy (|I8|I . The second condition would require that the minimum of the 
mapping curve is lower than the critical temperature of the original system without solvent (fig. |3J. 
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Figure 3: The dependence of the left hand side of (|14|l on the temperature for the MLG-type system 
(Mapping Curve) for 3 values of the protein-protein interaction energy. The horizontal dashed line shows 
the critical temperature of the original square well system (without solvent). All curves intersects the 
critical temperature of the system without solvent twice. This means that the MLG system in these 
cases has two critical points - one is the lower critical point (the slope of the curve is negative), the 
other is the upper critical point (the slope of the curve is positive). The phase diagram contains closed 
loop containing these two critical points. The parameters for the curves are the following E os = —2, 
E ob = -1, E d b = 1, E ds = 1.8, q os = 1, q ob = 10, q d b = 40, q ds = 49 



8 



4.3 The case in which the range of the solute-solute interactions, A, differs 
from the solute-solvent interactions, X s 

So far we have considered a special case in which the range of the solute-solute interaction is equal to 
the that of the shell region. In general this is not the case. In this section we develop an approximate 
mapping for such systems. We use Brilliantov's result 24 for the mean field relations between the critical 
temperature and critical density: 



= Pc 

v 



z + — 

2U4 



= h( Pc ) (28) 

CP 



where Vo 

v a = I v(r)dr (29) 



is the zeroth order moment of the attractive part v(r) of the potential. The terms inside the brackets on 
the right hand side of l)28[l are explained in |24| . This mean-field approximation has different precision 
for different interactions and even for the same interaction with different parameters. Division of the 
potential into the repulsive and attractive parts adds more freedom to obtaining parameter vq l|29l) , 
which also reduce quantitative precision. The assumption made in this section is that the critical density 
is the same for the systems with and without solvent. This is approximately true for A s < A = 1.25, where 
the critical density is about 0.4. For higher values of A this assumption fails. The analysis given in this 
section therefore can provide only qualitative understanding of the behavior of the critical temperature 
for different A s . 

What is important is that the right hand side of (|28|l depends only on the isothermal compressibility 
and its derivatives with respect to density of the hard sphere system PIJ. So for our purpose this is just 
some function of the density and independent of temperature. 

The square well system in the absence of solvent has vo equal to 

«o = -^(A 3 -l)e (30) 

where eo is the well depth of the particle particle interaction. For the system in the presence of the 
solvent i>o has the form: 

^ 3 /l3 ,x A^ms 1S _„ N , 2(e w - kT°As w ) Xl - 1 



(A - l)eo - 2(e w - kTAs w ) — ^—{X s -l) = «o - + =^ c ° - (31) 

o 3 \ e e A- 3 — I J 

Substituting vq and Vq into (|28|l and assuming that the critical density doesn't change we get the rela- 
tionship between the critical temperature of the system without solvent and the critical temperature of 
the system with solvent: 

"e , 2(e w -kT°As w )X 3 s -l~ 



fcTf = kT r 



(32) 



ksTcoex = — — — ; ^ 3 _i ; r(p c ) (33) 



A 3 - 1 

Solving i.Y3 for the kT^ we get the mapping relationship similar to (|15fl but for the critical temperature 
only: 

(eo + 2e w ) 
(l + 2As tu r( (0c )^i) 

For the special case A s equal to A this becomes l|15[l with r c = kT c /e. The main difference is that the 
relationship (|15fl is for all temperatures, while the l|33|l is only for the critical point. Let A° be the value 
of A s at which the denominator on the right hand side of l|33|) vanishes. For X s < A° the system has an 
upper critical point behavior. For X s > A° the system has a lower critical point behavior. 

Using (|32|) we can approximately map the phase diagram. In this case we can introduce the effective 
solvent parameters 

eff _ 

e ui — e w A 3 -l 



1 One of the methods to choose the form of the attractive part of potential is to set pair-correlation function g(r) equal 
to zero for r less than hard core diameter \'2h\ 
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These effective parameters correspond to the solute-soivent system with an effective shell size \ e Jf = A, 
with a zeroth moment vq equal to Vq. So we approximately change the system with solvent and A s =/= A 
by the system with solvent and effective parameters e%j* and As£/* with \ e /* = A. This allows us 
to perform a mapping of the SW system phase diagram onto the effective system with solvent phase 
diagram using (|32|l . We can consider the latter as the qualitative behavior of the system of interest. 

5 Modification of the Clausius-Clapeyron equation in the pres- 
ence of the solvent 

Now we consider a liquid-solid coexistence in the presence of the solvent. To determine the solid- liquid 
coexistence curve we use Gibbs-Duhem integration proposed by Kofke 2(i 27 . This method is based on 
the integration of the Clausius-Clapeyron equation: 

«f = Av (34) 

where AS and AV^ are the entropy and volume differences between solid and liquid phases, respectively. 
The derivative is taken along the coexistence curve. In the absence of the solvent this equation leads to 
the following form: 

dP _ Ae + PAv 

~dT ~ TAv {Ab) 
where Ae = e a — ei and Av are the differences in the energy and volume per particle in the solid and 
liquid phases respectively. 

When we have the model with a solvent contribution given by Q or (|5J), the interaction potential is 
temperature dependent. Therefore one must modify the (|35|) to take this temperature dependence into 
account. We start from the origins of the i|34|) . For the coexistent phases the difference of the Gibbs free 
energies is zero AG = G s — G; =0. Thus the derivative is taken at the conditions of AG = 



fdP\ _ (dAG/dT) NP _ AS 



(36) 



V^/ag=o (dAG/dP) NT AV 
ibs partition function into S = —dG/dT 

E{N, P, T) = C J dVdr N exp[-(3{U(r N , T) + PV)] (37) 



To get AS we substitute the Gibbs partition function into S = -dG/dT, where G = -k B T In S(JV, P, T), 
and 



The derivative of the Gibbs free energy with respect to temperature then is 
which leads to 

TS = -G+ < U > +P <V> -T < dU/dT > (38) 

Substituting this expression for the entropy into l|34(l we obtain the following modification of the Clausius- 
Clapeyron equation 

dP _ Ae + PAv - TAde 

dT ~ TAv ( } 

where de =< dU/dT > /N. 

When the Hamiltonian of the system consists of two body interactions between particles, the temper- 
ature is determined by the average kinetic energy and a temperature dependent microscopic Hamiltonian 
doesn't make much sense. In this case, however, we have a system of solute particles in a solvent environ- 
ment. A rearrangement of the solute particles leads to a rearrangement of solvent particles. The entropy 
of the system consists of two parts: the entropy of the protein particles and the entropy of the solvent 
molecules. Since in our model we don't consider separate solvent molecules, but rather just average their 
contribution, the effective interaction © has the temperature dependent term due to having integrated 
out the solvent degrees of freedom. This gives us the temperature dependent effective interaction. The 
last term in equation (|38|l can be considered as the solvent entropy contribution. 
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In the case of the square well protein-protein interactions, we can use the potential © . The derivative 
of this potential with respect to the temperature is 




01 t<><„. kBT dAs^_ kBAsw] , y V/) , (4()) 



dT \dT dT 

The —Nn c term in l|4l)|l cancels when we calculate the difference between solid and liquid phases. The 
average of this derivative per particle gives de = S(< n p > — n c ), where S is the expression in the first 
parenthesis in (|40|) and < n p > is the average number of protein-protein contacts. Using[5]we can relate 
the average number of contacts to the average energy. The difference between the average number of 
contacts in the solid and liquid phases is then related to the Ae: 



2Ae 

e + 2e w — 2k B TAs v 



A < n p >= — — — ~ (41) 



The difference of de between two phases is therefore 

A8e = 2SA : (42) 

e + 2e w - 2k B TAs w v ; 

Substituting ^ into (|2HJ> we obtain the final expression for the Clausius-Clapeyron equation for the 
square- well model with solvent 

d _P = / e + 2ew _ 2 ^ + 2kB T^ Ae \ i 

dT \ e + 2e w - 2k B TAs w / TAv v ; 



In the case of a system in which the solvent-solute interactions are not given by a square well potential 
we have 

de w i, rp dAs 

dT dT 



— = ( Ae + PAv 



k B T^ - k B A Sll 



e + 2e w — 2k B TAs v 



Ae-I^ (44) 



In equation l|35|) Ae, and Av are negative. The right hand side is therefore positive and the coexistent 
pressure increases with increasing temperature. In equation l|44|l Ae and At; are again negative, but the 
term in parenthesis can be either negative or positive; therefore for some cases the coexistent pressure can 
decrease as we increase the temperature. This effect is similar to the upside-down fluid-fluid coexistence 
curve and is due to the entropy contribution from the solvent. While the difference between the solid 
and the liquid of the entropy is generally negative (the solid state is more ordered) , the total solvent plus 
solute entropy difference can be positive. This means that the entropy increase of the solvent during 
protein crystallization due to the decrease of the contact number is greater than the entropy decrease of 
the solute. Therefore the numerator in i|34|) can be positive for some solvent parameters. At constant 
temperature AS* = AH/T. Thus this solvent model can reflect three different cases experimentally 
observed by Vekilov 0]. One is a negative enthalpy and entropy of crystallization, which gives the 
normal liquid-solid line. Another is a positive enthalpy of crystallization, which gives the upside down 
behavior of coexistence curves. The third is a zero enthalpy of crystallization, which gives the vertical 
coexistence curves. 



6 Numerical results for the solvent model 

6.1 Summary of the mapping procedure and results for the square well sys- 
tem 

In this section we summarize the algorithm of the mapping of the canonical square well (SW) phase 
diagram onto the diagram of the square well system in the presence of the solvent. The range of particle- 
particle interactions and the range of the particle-solvent interactions is considered to be the same. For 
the illustration we use the SW system with A = 1.25. Correspondingly X s = A = 1.25. The phase 
diagram of the system without solvent is shown in the fig. 21 The temperature of the system without 
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Figure 4: Phase diagram of the square well model with A = 1.25. Here 'x' denotes liquid-solidus 
line; '+' denotes fluid-fluid coexistence. The fluid- fluid coexistence line was obtained using the Gibbs 
ensemble method with N = 600 and V = 1500. The liquid-solid coexistence line was obtained using 
thermodynamic integration and Gibbs-Duhem integration techniques. Open square shows the position 
of the critical point. 



solvent we denote as r in order to not to confuse it with the temperature of the system with solvent, 
which is still denoted as kT. 

The dependence of the left hand side of (|14|) on the temperature of the system with solvent (kT in 
equation l|14|l 'l we will call the mapping curve. Now we construct the phase diagram of the system with 
the solvent. The mapping curve plays an important role in this construction process. Figure [S] shows the 
mapping procedure in detail. First we take the original phase diagram of the square well model(fig. [U<\). 
Second, choose some temperature n and invert the mapping curve (fig. |SJ3) to obtain the corresponding 
temperature of the system with the solvent 2 . Coexistent densities (pi and p4 in the fig. EJ) corresponding 
to the temperature t\ are also the coexistent densities of the system with solvent at a temperature kT\ . 
Now choose another temperature T2, invert it, using the mapping curve, to kT^ and use densities pi and 
P3 (fig- E3) as the coexistent densities at kT-i for the system with solvent. By continuing this process we 
obtain the phase diagram for the system with solvent (figure EJ). As one can see, this has a lower critical 
point and in general an upside down coexistence curve as compared with fig. 0] The interesting feature 
of this phase diagram is that at a temperature equal to 

»>— tst (45) 

the effective depth of the square well becomes zero and the particles behave as hard spheres, with the 
corresponding coexistent densities equal to the hard spheres case. In the figure [5] kTjj 5 = 1/3. 

The next step is to invert the mapping curves shown in the figure [21 For this purpose we have to 
solve numerically equation (|14[) with parameters given by ill'Ot - l|23l) . Figure shows three fluid-fluid 
coexistence curves corresponding to the three mapping curves shown on fig. One can see that indeed 
the phase diagrams have the form of closed loops. This is because each temperature r of the system 
without solvent corresponds to two temperatures kT\ and kTi for the system with solvent. This can be 
understood from consideration of the MLG energy level model (fig. [5J . At low temperatures the solvent 

2 Note that for constant e w and As w this inversion has an analytical solution 1151 . However for temperature dependent 
e w and As w as in the case of the closed loop like phase diagram, we have to invert the mapping curve numerically. 
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Figure 5: Mapping scheme. A. The fluid-fluid coexistence line for the square well model in the absence 
of solvent (same as on fig. 0J. B. The mapping curve for e w = —1 and As w = —1.5. C. Results of the 
mapping. Filled squares are the points on the corresponding phase diagram for the square well model 
with solvent. See text for more description. 
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Figure 6: The phase diagram of the square well system with the solvent parameters e w = — 1 and 
Asu, = —1.5, obtained by the procedure shown in figureEl 'x' denotes liquid-solidus line; '+' denotes 
fluid-fluid coexistence. 
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Figure 7: Fluid-fluid coexistence curves in the form of closed loop for the different protein-protein 
interaction strength e . The choice of solvent parameters is determined by equations <{^0|) - The 
case with eo = qualitatively corresponds to the situation described in pQ where two proteins interact 
as hard spheres (on lattice). 



molecules mostly occupy the lowest level, corresponding to the ordered shell. This state corresponds to 
the uniform distribution of the solute particles surrounded by the structures of solvent molecules. At 
high temperatures the solvent molecules mostly occupy the top energy level, which corresponds to the 
disordered shell. So the system again tends to be in the uniform state, but now the solvent surrounding 
the solute particles is disordered. At intermediate temperature the system tends to occupy the bulk 
energy levels, which favors the phase separation. This is the property of the MLG model that reflects 
the picture described in the introduction. 

The mapping curves shown in the figure [21 produce the closed loop like fluid-fluid coexistence curves. 
Figure [H] shows the behavior of the liquid-solid line in this case. The liquid-solid line also displays a 
behavior similar to the fluid fluid coexistence line. We are unaware of any experimental observations of 
such a liquid-solid coexistence behavior and do not know whether it is simply an artifact of our particular 
model. 

6.2 Numerical results for the MLJ system in the presence of solvent 

In this section consider the effect of solvent on particles interacting via a modified Lennard- Jones (MLJ) 
potential given by 

V(r) = { 4e f i ° °' a ) [ > I (46) 

It has been shown that at the critical point the nucleation rate is many orders of magnitude greater 
than at other points in the phase diagram. This suggests that the nucleation of protein particles can 
be achieved near the critical point. The model has been well studied and its critical point accurately 
determined. 

For the case of the MLJ model, we again choose to use a solvent-solute interaction range, X s , that is 
equal to the protein-protein particle interaction range. Because the particles interacting via this potential 
have an effective hard-core diameter, care must be taken such that we choose an appropriate value for 



14 



1.6 



1.4 



0.6 





X 


X 




X 


X 




X 


X 


- 


X 


X 




X 


X 
X 




X 


X 




X 


X 




X 


X 


X 


X 


X 
X 


- /> ++ + 

x rf*- 
x t 
x + 

x + 


+ 

+ 


X 
X 

X 


x + 

w 

x+ 


+ 


X 






x 

X *x x 

X 

X 

>k 






X 
X 



0.2 0.4 0.6 0.8 1 1.2 

P 

Figure 8: The fluid- fluid coexistence line as well as the solid-liquid coexistence line for the protein-protein 
interaction strength eo — for the choice of parameters determined by equations l|20|) - (|25|l . 



A s . We calculate the effective hard-core diameter using 

/>oo 

<7e//= / dr[l-exp{-Vrep/k B T)], (47) 
Jo 

where V rep is the repulsive part of the potential. The range over which the particles interact is controlled 
by the parameter a in 16.21 We choose a = 50 as in other studies. It has been shown that for a = 50, 
one can obtain an equivalent range in terms of A, the parameter denoting the range in the square well 
system. We use the value A = 1.073. To account for the effects of the effective hard-core diameter, we 
actually use A s = 1.26 for the solvent-solute interaction range. The values e w = —1 and AS W = —1.5 
were used as before to obtain an upside-down phase diagram. 

To calculate the solid-fluid phase boundaries, we use a modified Gibbs-Duhem equation given by eq. 
43. A coexistence point was calculated using free-energy methods and simulations of a coupled Einstein 
lattice. Isobaric-isothermal (NPT) simulations were performed in parallel for N = 256 particles on a 
periodic simulation cell to obtain the entire coexistence curve. Equilibration and production times were 
five million and ten million Monte Carlo steps, respectively. 

We calculated the fluid-fluid coexistence curve using the Gibbs ensemble Monte Carlo method. Two 
physically separated, but thermodynamically connected, simulation cells are allowed to exchange particles 
and undergo volume displacements such that the total number of particles N = N1+N2 and total volume 
V = V1 + V2 remain constant. Simulations were performed on an N = 600 particle system. Equilibration 
and production runs were fifty million and one hundred million Monte Carlo steps, respectively. Our 
results are shown in figure [§J 

As can be seen from the figure, the phase diagram for this model with the particular parameters used 
is very similar to that for the square well model shown in figure [6] 



7 Conclusion 

The model presented in this paper takes into account the solvent contribution to the solute-solvent free 
energy of globular proteins in solution. The contribution depends on the parameters that describe the 
free energy of the solvent molecule change, e w and As w . These parameters play the role of the solvent 
enthalpy and entropy change per solvent particle as the molecule goes from the bulk to the vicinity of 
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Figure 9: Phase diagram of the ML J model including solvent-solute interactions. Interactions between the latter 
were mediated over the range A s = 1.26. Again 'x' denotes liquid-solidus line; '+' denotes fluid-fluid coexistence. 



the protein molecule. The solvent enthalpy and entropy change per solute particle upon changing from 
the fluid phase to the solid phase(parameters AH so i ven t and AS so ivent from can be calculated as 
equal to e w and As w times the difference of the average number of contacts per molecule in these two 
phases. 

AH solvent = e w (< 4 lmd >-< n s ° Hte >) 
AS solvent = As w (< n{ luld > - < n s c ° lute >) 

These relationships allow us to relate e w and As w to the three cases of the solubility dependence 
described in the introduction - normal solubility dependence, retrograde solubility dependence and con- 
stant solubility (as in the case of apoferritin) . All we need is the sign of the total enthalpy change AH. 
For the square well case, e w > — eo/2 is a condition for AH being negative and therefore the solubility 
curve being normal. If the range of attraction is short this qualitatively describes the lysozyme phase 
diagram. The condition e w < —eo/2 is a condition for AH being positive and therefore the solubility 
curve being retrograde. This corresponds to the case of the HbC solubility curve (fig. Q). Indeed we can 
see that the liquidus line in the figure[S]has qualitatively the same behavior as the HbC solubility curve. 
If e w = e/2, so that the enthalpy change is small, then the solubility curve is almost vertical, with the 
sign of the slope determined by the sign of the enthalpy change. 

The advantage of this simplified model is that in the particular case of the square well, with the 
particle-particle range of interaction equal to the width of the shell region around the particle, one can 
obtain the phase diagram by a mapping of the square well phase diagram without solvent. Therefore 
one doesn't have to perform Monte Carlo simulations or theoretical approximations to obtain the phase 
diagram for the model with solvent. 

The exact mapping is also possible for the hard sphere model with solvent. In this case eo = 0. The 
hard sphere case is interesting in two aspects. First, one can obtain an exact mapping procedure for any 
solvent-solute interaction. Second, the hard sphere interaction doesn't include any initial solvent contri- 
bution (unlike the case with an effective solute-solute interaction). This can be the case of noninteracting 
neutral colloidal particles in the solvent. If we put such particles into water, the hydrogen bonds break 
and rearrange, as explained in the introduction and in pp. Therefore we can use the MLG model to get 
solvent parameters e w and As w . We can consider the hard sphere interaction to be a square well with 
eo = and A = X s . Then using a mapping curve as in the figure 01 (the one with eo = 0) we can obtain 
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the phase diagram of the neutral, noninteracting spherical particles in water. If the critical temperature 
of the square well system with a range of interaction equal to A s intersects the mapping curve (fig. |3J), the 
colloidal system has a closed loop type phase diagram (fig with eo = 0). As the value of A s decreases 
(the size of the particle increases), the fluid-fluid coexistence curve shrinks and disappears. So larger 
noninteracting particles don't have a fluid-fluid phase separation in this model. However, if the particles 
starts to interact, the mapping curve lowers (fig[3J) and phase separation may occur again. 

For the case of constant solvent parameters we can have either an upper critical point with a normal 
solubility dependence, or a lower critical point with a retrograde solubility dependence. The temperature 
dependent parameters allow us to have a combination of these behaviors. A natural way to derive the 
temperature dependent solvent parameters for water is to use the MLG model. This model, together 
with our solvent-solute interaction, has closed loop phase diagrams, similar to £Q. 
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